R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

library(xgboost)

载入程辑包：‘xgboost’

The following object is masked from ‘package:IRanges’:

    slice

The following object is masked from ‘package:plotly’:

    slice

The following object is masked from ‘package:dplyr’:

    slice

library(Matrix)

载入程辑包：‘Matrix’

The following object is masked from ‘package:S4Vectors’:

    expand

The following objects are masked from ‘package:tidyr’:

    expand, pack, unpack

library(mclust)

    __  ___________    __  _____________
   /  |/  / ____/ /   / / / / ___/_  __/
  / /|_/ / /   / /   / / / /\__ \ / /   
 / /  / / /___/ /___/ /_/ /___/ // /    
/_/  /_/\____/_____/\____//____//_/    version 5.4.9
Type 'citation("mclust")' for citing this R package in publications.

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
─ Attaching packages ───────────────────────────────────────── tidyverse 1.3.1 ─
✓ tibble  3.1.5     ✓ stringr 1.4.0
✓ readr   2.0.2     ✓ forcats 0.5.1
✓ purrr   0.3.4     
─ Conflicts ────────────────────────────────────────── tidyverse_conflicts() ─
x IRanges::collapse()      masks dplyr::collapse()
x Biobase::combine()       masks BiocGenerics::combine(), dplyr::combine()
x IRanges::desc()          masks dplyr::desc()
x Matrix::expand()         masks S4Vectors::expand(), tidyr::expand()
x plotly::filter()         masks dplyr::filter(), stats::filter()
x S4Vectors::first()       masks dplyr::first()
x widgetTools::funs()      masks dplyr::funs()
x dplyr::lag()             masks stats::lag()
x purrr::map()             masks mclust::map()
x Matrix::pack()           masks tidyr::pack()
x BiocGenerics::Position() masks ggplot2::Position(), base::Position()
x purrr::reduce()          masks IRanges::reduce()
x S4Vectors::rename()      masks plotly::rename(), dplyr::rename()
x AnnotationDbi::select()  masks plotly::select(), dplyr::select()
x purrr::simplify()        masks clusterProfiler::simplify()
x xgboost::slice()         masks IRanges::slice(), plotly::slice(), dplyr::slice()
x Matrix::unpack()         masks tidyr::unpack()

ds0 <- readRDS("./ds0.rds")
ds1 <- readRDS("./ds1.rds")

Loading required package: Seurat
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
Registered S3 method overwritten by 'spatstat':
  method     from
  print.boxx cli 

ds2 <- readRDS("./ds2.rds")

分发训练集

# Idents(ds2) <- ds2$conditions
# ds2_AC <- subset(ds2, idents = "AC")
# ds2_PA <- subset(ds2, idents = "PA")
# ds2_AC <- ds2_AC %>% FindNeighbors(dims = 1:20) %>% FindClusters(resolution = 0.1)
# ds2_PA <- ds2_PA %>% FindNeighbors(dims = 1:20) %>% FindClusters(resolution = 0.1)
# umapplot(ds2_AC) + scale_y_continuous(limits = c(-5,15),breaks = NULL) +
#         scale_x_continuous(limits = c(-5,15),breaks = NULL)
# umapplot(ds2_PA)+ scale_y_continuous(limits = c(-5,15),breaks = NULL) +
#         scale_x_continuous(limits = c(-5,15),breaks = NULL)
# 
# AC_markers <- FindAllMarkers(ds2_AC,logfc.threshold = 0.7, min.diff.pct = 0.2)
# # PA_markers <- FindAllMarkers(ds2_PA,logfc.threshold = 0.7, min.diff.pct = 0.2)
# write.csv(AC_markers,"AC_SMC_markers.csv")

ds2_AC <- readRDS("ds2_AC.rds")
ds2_PA <- readRDS("ds2_PA.rds")

umapplot(ds2_AC)
umapplot(ds2_PA)

ds2_AC$Classification <- Idents(ds2_AC)
Idents(ds2_AC) <- ds2_AC$seurat_clusters
ds2_AC <- RenameIdents(ds2_AC,'0' = '3','1' = '1','2' = '0','3' = '2')
Idents(ds2_AC) <- factor(Idents(ds2_AC),levels = c(0,1,2,3))
ds2_AC$seurat_clusters <- Idents(ds2_AC)

ds2_PA$Classification <- Idents(ds2_PA)
Idents(ds2_PA) <- ds2_PA$seurat_clusters

在AC上预训练

ds2_AC$Classification <- Idents(ds2_AC)
Idents(ds2_AC) <- ds2_AC$seurat_clusters
AC_data <- get_data_table(ds2_AC, highvar = F, type = "data")
AC_label <- as.numeric(as.character(Idents(ds2_AC)))

set.seed(7)
index <- c(1:dim(AC_data)[2]) %>% sample(ceiling(0.3*dim(AC_data)[2]), replace = F, prob = NULL)

colnames(AC_data) <- NULL

AC_train_data <- list(data = t(as(AC_data[,-index],"dgCMatrix")), label = AC_label[-index])
AC_test_data <- list(data = t(as(AC_data[,index],"dgCMatrix")), label = AC_label[index])

AC_train <- xgb.DMatrix(data = AC_train_data$data,label = AC_train_data$label)
AC_test <- xgb.DMatrix(data = AC_test_data$data,label = AC_test_data$label)

# xgb_params_train = {
#     'objective':'multi:softmax',
#     'eval_metric':'mlogloss',
#     'num_class':self.numbertrainclasses,
#     'eta':0.2,
#     'max_depth':6,
#     'subsample': 0.6}
# nround = 200

watchlist <- list(train = AC_train, eval = AC_test)
xgb_param <- list(eta = 0.2, max_depth = 6, 
                  subsample = 0.6,  num_class = length(table(Idents(ds2_AC))),
                  objective = "multi:softmax", eval_metric = 'mlogloss')

bst_model <- xgb.train(xgb_param, AC_train, nrounds = 100, watchlist, verbose = 0)

eval_loss <- bst_model[["evaluation_log"]][["eval_mlogloss"]]
plot_ly(data.frame(eval_loss), x = c(1:100), y = eval_loss) %>% 
  add_trace(type = "scatter", mode = "markers+lines", 
            marker = list(color = "black", line = list(color = "#1E90FFC7", width = 1)),
            line = list(color = "#1E90FF80", width = 2)) %>% 
  layout(xaxis = list(title = "epoch"),yaxis = list(title = "eval_mlogloss"))

# 特征提取
importance <- xgb.importance(colnames(AC_train), model = bst_model)
head(importance)
xgb.ggplot.importance(head(importance,20), n_clusters = 1) + theme_minimal()

multi_featureplot(head(importance,9)$Feature, ds2_AC) 
AC_genes <- head(importance, 500) ##选择top500

write.csv(AC_genes, "./datatable/AC_features.csv", row.names = F)

#混淆矩阵
predict_AC_test <- round(predict(bst_model, newdata = AC_test))

AC_confuse_matrix_test <- table(AC_test_data$label, predict_AC_test, dnn=c("true","pre"))
AC_confuse_matrix_test_prop <- prop.table(AC_confuse_matrix_test, 1)
AC_confuse_matrix_test_prop

confuse_bubblemat(AC_confuse_matrix_test_prop, c("Fibroblast", "SMC1", "Fibromyocyte", "SMC2"), c("Fibroblast", "SMC1", "Fibromyocyte", "SMC2"),"AC_pretrain")

#ROC曲线
xgboost_roc <- pROC::multiclass.roc(AC_test_data$label, predict_AC_test) #多分类ROC
xgboost_roc[["auc"]] #只需要这个值
adjustedRandIndex(AC_test_data$label, predict_AC_test) #分类器性能

在PA上训练

ds2_PA$Classification <- Idents(ds2_PA)
Idents(ds2_PA) <- ds2_PA$seurat_clusters

PA_data <- get_data_table(ds2_PA, highvar = F, type = "data")
PA_label <- as.numeric(as.character(Idents(ds2_PA)))
set.seed(7)
index <- c(1:dim(PA_data)[2]) %>% sample(ceiling(0.3*dim(PA_data)[2]), replace = F, prob = NULL)
colnames(PA_data) <- NULL

PA_train_data <- list(data = t(as(PA_data[,-index],"dgCMatrix")), label = PA_label[-index])
PA_test_data <- list(data = t(as(PA_data[,index],"dgCMatrix")), label = PA_label[index])

PA_train <- xgb.DMatrix(data = PA_train_data$data,label = PA_train_data$label)
PA_test <- xgb.DMatrix(data = PA_test_data$data,label = PA_test_data$label)

watchlist <- list(train = PA_train, eval = PA_test)
xgb_param <- list(eta = 0.2, max_depth = 6, 
                  subsample = 0.6,  num_class = length(table(Idents(ds2_PA))),
                  objective = "multi:softmax", eval_metric = 'mlogloss')
bst_model <- xgb.train(xgb_param, PA_train, nrounds = 100, watchlist, verbose = 0)
eval_loss <- bst_model[["evaluation_log"]][["eval_mlogloss"]]
plot_ly(data.frame(eval_loss), x = c(1:100), y = eval_loss) %>% 
  add_trace(type = "scatter", mode = "markers+lines", 
            marker = list(color = "black", line = list(color = "#1E90FFC7", width = 1)),
            line = list(color = "#1E90FF80", width = 2)) %>% 
  layout(xaxis = list(title = "epoch"),yaxis = list(title = "eval_mlogloss"))

# 特征提取
importance <- xgb.importance(colnames(PA_train), model = bst_model)
head(importance)
xgb.ggplot.importance(head(importance,20),n_clusters = 1) + theme_minimal()

multi_featureplot(head(importance,9)$Feature, ds2_PA)
PA_genes <- head(importance, 500) ##选择top500
write.csv(PA_genes, "./datatable/PA_features.csv", row.names = F)

#混淆矩阵
predict_PA_test <- round(predict(bst_model, newdata = PA_test))

PA_confuse_matrix_test <- table(PA_test_data$label, predict_PA_test, dnn=c("true","pre"))
PA_confuse_matrix_test_prop <- prop.table(PA_confuse_matrix_test,1)
PA_confuse_matrix_test_prop

confuse_bubblemat(PA_confuse_matrix_test_prop,c("Fibromyocyte", "SMC1", "SMC2"),c("Fibromyocyte", "SMC1", "SMC2"),"PA_pretrain")

#ROC曲线

xgboost_roc <- pROC::multiclass.roc(PA_test_data$label, predict_PA_test) #多分类ROC
xgboost_roc[["auc"]]
adjustedRandIndex(PA_test_data$label, predict_PA_test) #PA分类器性能

选择特征common genes of top 500

使用所有来自PA的细胞训练分类器

应用在AC上，计算ARI

selected_features <- read.csv("./datatable/selected_features.csv", stringsAsFactors = F)
selected_features <- selected_features$x
PA_data <- get_data_table(ds2_PA, highvar = F, type = "data")
PA_data <- PA_data[selected_features,]
PA_label <- as.numeric(as.character(Idents(ds2_PA)))
colnames(PA_data) <- NULL

PA_train_data <- list(data = t(as(PA_data,"dgCMatrix")), label = PA_label)
PA_train <- xgb.DMatrix(data = PA_train_data$data,label = PA_train_data$label)
xgb_param <- list(eta = 0.2, max_depth = 6, 
                  subsample = 0.6,  num_class = length(table(Idents(ds2_PA))),
                  objective = "multi:softmax", eval_metric = 'mlogloss')

bst_model <- xgb.train(xgb_param, PA_train, nrounds = 100, verbose = 0)

# 特征提取
importance <- xgb.importance(colnames(PA_train), model = bst_model)
head(importance)

xgb.ggplot.importance(head(importance,20),n_clusters = 1) + theme_bw()

write.csv(importance, "./datatable/PAtrain_features.csv", row.names = F)

应用到AC上

AC_data <- get_data_table(ds2_AC, highvar = F, type = "data")
AC_data <- AC_data[selected_features,]
AC_label <- as.numeric(as.character(Idents(ds2_AC)))
colnames(AC_data) <- NULL
AC_test_data <- list(data = t(as(AC_data,"dgCMatrix")), label = AC_label)
AC_test <- xgb.DMatrix(data = AC_test_data$data,label = AC_test_data$label)

#计算混淆矩阵
predict_AC_test <- round(predict(bst_model, newdata = AC_test))
AC_confuse_matrix_test <- table(AC_test_data$label, predict_AC_test, dnn=c("true","pre"))
AC_confuse_matrix_test_prop <- prop.table(AC_confuse_matrix_test,1)
AC_confuse_matrix_test_prop  #分析发育轨迹

    pre
true           0           1           2
   0 0.996726678 0.003273322 0.000000000
   1 0.825876663 0.170495768 0.003627570
   2 0.435185185 0.509259259 0.055555556
   3 0.002762431 0.069060773 0.928176796

confuse_bubblemat(AC_confuse_matrix_test_prop,c("Fibroblast", "SMC1", "Fibromyocyte", "SMC2"), c("Fibromyocyte", "SMC1", "SMC2"), "PAtoAC")


#ROC曲线
xgboost_roc <- pROC::multiclass.roc(AC_test_data$label, predict_AC_test) #多分类ROC

Setting direction: controls < cases
Setting direction: controls < cases
Setting direction: controls < cases
Setting direction: controls < cases
Setting direction: controls < cases
Setting direction: controls < cases

xgboost_roc[["auc"]]

Multi-class area under the curve: 0.8342

# 计算ARI 
adjustedRandIndex(predict_AC_test, AC_test_data$label)

[1] 0.3147662

sankey plot

PA -> AC

sankey_plot(AC_confuse_matrix_test, label1 = c("Fibroblast", "SMC1", "SMC2"), label2 = c("Fibroblast", "SMC1", "Fibromyocyte", "SMC2"), session = "PA -> AC")

umapplot(ds2_AC)

umapplot(ds2_PA)

# umapplot(ds2,split.by = "conditions")

#把结果投射回umap

Idents(ds2_AC) <- predict_AC_test
ds2_AC$predict_AC_test <- predict_AC_test
umapplot(ds2_AC,group.by = "predict_AC_test")

Idents(ds2_AC) <- ds2_AC$seurat_clusters

反着做

选择特征common genes of top 500

使用所有来自AC的细胞训练分类器

应用在PA上，计算ARI

PA_data <- get_data_table(ds2_PA, highvar = F, type = "data")
PA_data <- PA_data[selected_features,]
PA_label <- as.numeric(as.character(Idents(ds2_PA)))
colnames(PA_data) <- NULL

PA_test_data <- list(data = t(as(PA_data,"dgCMatrix")), label = PA_label)

PA_test <- xgb.DMatrix(data = PA_test_data$data,label = PA_test_data$label)

#计算混淆矩阵
predict_PA_test <- round(predict(bst_model2, newdata = PA_test))
 
PA_confuse_matrix_test <- table(PA_test_data$label, predict_PA_test, dnn=c("true","pre"))
PA_confuse_matrix_test_prop <- prop.table(PA_confuse_matrix_test,1)
PA_confuse_matrix_test_prop  #分析发育轨迹

confuse_bubblemat(PA_confuse_matrix_test_prop,c("Fibromyocyte", "SMC1", "SMC2"),c("Fibroblast", "SMC1", "Fibromyocyte", "SMC2"),session = "ACtoPA")

# 计算ARI
adjustedRandIndex(predict_PA_test, PA_test_data$label)

把结果投射回umap

Idents(ds2_PA) <- predict_PA_test
ds2_PA$predict_PA_test <- predict_PA_test
umapplot(ds2_PA,group.by = "predict_PA_test")
Idents(ds2_PA) <- ds2_PA$seurat_clusters

sankey plot

labels <- c("Fibroblast", "SMC1", "Fibromyocyte", "SMC2")
labels2 <- c("Fibromyocyte", "SMC1", "SMC2")
sankey_plot(PA_confuse_matrix_test,labels,labels2,session = "AC -> PA")

varify 部分

病变程度量化 # 数据集CA_dataset1 ## ds2全体训练

Idents(ds2) <- ds2$seurat_clusters
ds2_data <- get_data_table(ds2, highvar = F, type = "data")
ds2_label <- as.numeric(as.character(Idents(ds2)))

index <- c(1:dim(ds2_data)[2]) %>% sample(ceiling(0.3*dim(ds2_data)[2]), replace = F, prob = NULL)
colnames(ds2_data) <- NULL

ds2_train_data <- list(data = t(as(ds2_data[,-index],"dgCMatrix")), label = ds2_label[-index])
ds2_test_data <- list(data = t(as(ds2_data[,index],"dgCMatrix")), label = ds2_label[index])

ds2_train <- xgb.DMatrix(data = ds2_train_data$data,label = ds2_train_data$label)
ds2_test <- xgb.DMatrix(data = ds2_test_data$data,label = ds2_test_data$label)

watchlist <- list(train = ds2_train, eval = ds2_test)
xgb_param <- list(eta = 0.2, max_depth = 6, 
                  subsample = 0.6,  num_class = length(table(Idents(ds2))),
                  objective = "multi:softmax", eval_metric = 'mlogloss')

bst_model <- xgb.train(xgb_param, ds2_train, nrounds = 100, watchlist, verbose = 0)

eval_loss <- bst_model[["evaluation_log"]][["eval_mlogloss"]]
plot_ly(data.frame(eval_loss), x = c(1:100), y = eval_loss) %>% 
  add_trace(type = "scatter", mode = "markers+lines", 
            marker = list(color = "black", line = list(color = "#1E90FFC7", width = 1)),
            line = list(color = "#1E90FF80", width = 2)) %>% 
  layout(xaxis = list(title = "epoch"),yaxis = list(title = "eval_mlogloss"))

# 特征提取
importance <- xgb.importance(colnames(ds2_train), model = bst_model)
head(importance)
xgb.ggplot.importance(head(importance,20), n_clusters = 1) + theme_minimal()

multi_featureplot(head(importance,9)$Feature, ds2) 
ds2_genes <- head(importance, 500) ##选择top500

write.csv(ds2_genes, "./datatable/ds2_features.csv", row.names = F)


predict_ds2_test <- round(predict(bst_model, newdata = ds2_test)) %>% 
#混淆矩阵
ds2_confuse_matrix_test <- table(ds2_test_data$label, predict_ds2_test, dnn=c("true","pre"))
ds2_confuse_matrix_test_prop <- prop.table(ds2_confuse_matrix_test, 1)
ds2_confuse_matrix_test_prop

x <- c(0:4)
y <- c(0:4)
confuse_bubblemat(ds2_confuse_matrix_test_prop,x,y,"ds2_train")


#ROC曲线
xgboost_roc <- pROC::multiclass.roc(ds2_test_data$label, predict_ds2_test) #多分类ROC
xgboost_roc[["auc"]] #只需要这个值
adjustedRandIndex(ds2_test_data$label, predict_ds2_test) #分类器性能

temp <- get_data_table(ds1, highvar = T, type = "data")
ds1_data <- matrix(data=0, nrow = dim(ds2_data)[1], ncol = length(colnames(temp)), 
                   byrow = FALSE, dimnames = list(rownames(ds2_data),colnames(temp)))
for(i in intersect(rownames(ds2_data), rownames(temp))){
  ds1_data[i,] <- temp[i,]
}
rm(temp)
ds1_label <- as.numeric(as.character(Idents(ds1)))
colnames(ds1_data) <- NULL
ds1_test_data <- list(data = t(as(ds1_data,"dgCMatrix")), label = ds1_label)
ds1_test <- xgb.DMatrix(data = ds1_test_data$data,label = ds1_test_data$label)

#预测结果
predict_ds1_test <- round(predict(bst_model, newdata = ds1_test))

#计算混淆矩阵
ds1_data_confuse_matrix_test <- table(ds1_test_data$label, predict_ds1_test, dnn=c("true","pre"))
ds1_data_confuse_matrix_test_prop <- prop.table(ds1_data_confuse_matrix_test,1)

#绘制混淆矩阵
x <- c("Fibromyocyte", "SMC1", "SMC2")
y <- c("Fibroblast", "SMC1", "Fibromyocyte", "SMC2")
confuse_bubblemat(ds1_data_confuse_matrix_test_prop,x,y,"ds2tods1")

ds1_data_confuse_matrix_test
ds1_data_confuse_matrix_test_prop  #分析发育轨迹
#ROC曲线
xgboost_roc <- pROC::multiclass.roc(ds1_test_data$label, predict_ds1_test) #多分类ROC
xgboost_roc[["auc"]]

# 计算ARI 
adjustedRandIndex(predict_ds1_test, ds1_test_data$label)

投射回umap

Idents(ds1) <- predict_ds1_test
ds1$predict_ds1_test <- predict_ds1_test
umapplot(ds1,group.by = "predict_ds1_test")
Idents(ds1) <- ds1$seurat_clusters

冠状动脉数据集

ds0 <- ds0 %>% FindNeighbors(dims = 1:20) %>% FindClusters(resolution = 0.1)
umapplot(ds0)
f("TAGLN",ds0)
# ds0_markers <- FindAllMarkers(ds0,logfc.threshold = 0.7, min.diff.pct = 0.2)

selected_features <- AC_genes$Feature
temp <- get_data_table(ds0, highvar = F, type = "data")
ds0_data <- matrix(data=0, nrow = length(selected_features), 
                   ncol = length(colnames(temp)), byrow = FALSE, 
                   dimnames = list(selected_features,colnames(temp)))
for(i in intersect(selected_features,rownames(temp))){
  ds0_data[i,] <- temp[i,]
}
rm(temp)

ds0_label <- as.numeric(as.character(Idents(ds0)))
colnames(ds0_data) <- NULL
ds0_test_data <- list(data = t(as(ds0_data,"dgCMatrix")), label = ds0_label)
ds0_test <- xgb.DMatrix(data = ds0_test_data$data,label = ds0_test_data$label)

#计算混淆矩阵
predict_ds0_test <- round(predict(bst_model, newdata = ds0_test))

ds0_data_confuse_matrix_test <- table(ds0_test_data$label, predict_ds0_test, dnn=c("true","pre"))
ds0_data_confuse_matrix_test_prop <- prop.table(ds0_data_confuse_matrix_test,1)

x <- c("ds0_0", "ds0_1", "ds0_2", "ds0_3", "ds0_4")
y <- c("AC_0", "AC_1", "AC_2")

prop <- as.numeric(ds0_data_confuse_matrix_test_prop)
data <- expand.grid(x = x, y = y) %>% bind_cols(prop = prop)
plot <- ggplot(data, aes(x = x, y = y, colour = prop, size = prop)) +
  geom_point()+
  scale_size_continuous(range = c(0, 10)) + 
  labs(x = "clusters", y = "inferred from") + theme_bw()
ggsave("./plots/ACmodel_humancor.png", plot = plot, device = png, width = 5,height = 4)

ds0_data_confuse_matrix_test
ds0_data_confuse_matrix_test_prop  #分析发育轨迹

#ROC曲线
xgboost_roc <- pROC::multiclass.roc(ds0_test_data$label, predict_ds0_test) #多分类ROC

# 计算ARI 
adjustedRandIndex(predict_ds0_test, ds0_test_data$label)

labels <- lapply(levels(Idents(ds2_AC)), paste0, "_AC") %>% as.character()
labels2 <- lapply(levels(Idents(ds0)), paste0, "_ds0") %>% as.character()
sources <- rep(0:(length(labels)-1), each = length(labels2))  #注意这里的each和times的区别
colors <- rep(colors_list[1:length(labels)], each = length(labels2))
targets <- rep(length(labels)+0:(length(labels2)-1), times = length(labels))

plot_ly(type = "sankey", orientation = "h",
    node = list(
      label = c(labels,labels2), 
      color = colors_list, pad = 15, thickness = 30,
      line = list(
        color = "black",
        width = 1)),
    link = list(
      source = sources,
      target = targets,
      value =  as.numeric(ds0_data_confuse_matrix_test),
      color = colors
      ))

# load("./init.RData")
multi_featureplot(head(importance2,9)$Feature, ds2_AC)
multi_featureplot(head(importance2,9)$Feature, ds0)
multi_featureplot(head(importance2,9)$Feature, ds1)
f("MYH11", ds2_AC)
umapplot(ds0)

淋巴细胞

Idents(lym_ds2) <- lym_ds2$conditions
lym_ds2_AC <- subset(lym_ds2, idents = "AC")
lym_ds2_PA <- subset(lym_ds2, idents = "PA")
lym_ds2_AC <- lym_ds2_AC %>% FindNeighbors(dims = 1:20) %>% FindClusters(resolution = 0.2)

Computing nearest neighbor graph
Computing SNN

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 2990
Number of edges: 98189

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.9128
Number of communities: 5
Elapsed time: 0 seconds

umapplot(lym_ds2_AC)

lym_ds2_PA <- lym_ds2_PA %>% FindNeighbors(dims = 1:20) %>% FindClusters(resolution = 0.2)

Computing nearest neighbor graph
Computing SNN

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 13746
Number of edges: 456548

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.9116
Number of communities: 6
Elapsed time: 2 seconds

umapplot(lym_ds2_PA)

用PA的lym训练

lym_PA_data <- get_data_table(lym_ds2_PA, highvar = F, type = "data")
lym_PA_label <- as.numeric(as.character(Idents(lym_ds2_PA)))

set.seed(7)
index <- c(1:dim(lym_PA_data)[2]) %>% sample(ceiling(0.3*dim(lym_PA_data)[2]), replace = F, prob = NULL)
colnames(lym_PA_data) <- NULL
lym_PA_train_data <- list(data = t(as(lym_PA_data[,-index],"dgCMatrix")), label = lym_PA_label[-index])
lym_PA_test_data <- list(data = t(as(lym_PA_data[,index],"dgCMatrix")), label = lym_PA_label[index])

lym_PA_train <- xgb.DMatrix(data = lym_PA_train_data$data,label = lym_PA_train_data$label)
lym_PA_test <- xgb.DMatrix(data = lym_PA_test_data$data,label = lym_PA_test_data$label)

watchlist <- list(train = lym_PA_train, eval = lym_PA_test)
xgb_param <- list(eta = 0.2, max_depth = 6, 
                  subsample = 0.6,  num_class = length(table(Idents(lym_ds2_PA))),
                  objective = "multi:softmax", eval_metric = 'mlogloss')

bst_model <- xgb.train(xgb_param, lym_PA_train, nrounds = 100, watchlist, verbose = 0)

用AC的lym验证

labels <- lapply(levels(Idents(lym_ds2_PA)), paste0, "_lymPA") %>% as.character()
labels2 <- lapply(levels(Idents(lym_ds2_AC)), paste0, "_lymAC") %>% as.character()
sources <- rep(0:5, each = 5)  #注意这里的each和times的区别
colors <- rep(colors_list[1:6], each = 5)
targets <- rep(6:10, times = 6)

plot_ly(type = "sankey", orientation = "h",
    node = list(
      label = c(labels,labels2), 
      color = colors_list, pad = 15, thickness = 30,
      line = list(
        color = "black",
        width = 1)),
    link = list(
      source = sources,
      target = targets,
      value =  as.numeric(lym_AC_confuse_matrix_test),
      color = colors
      ))


umapplot(lym_ds2_AC)
umapplot(lym_ds2_PA)

functions set

sankey_plot <- function(confuse_matrix, label1, label2, session = "session")
{
  sources <- rep(0:(length(label1)-1), each = length(label2))  #注意这里的each和times的区别
  colors <- rep(aero_colors_list[1:length(label1)], each = length(label2))
  targets <- rep(length(label1)+0:(length(label2)-1), times = length(label1))

  plot_ly(type = "sankey", orientation = "h",
      node = list(
        label = c(label1,label2), 
        color = colors_list, pad = 15, thickness = 30,
        line = list(color = "black", width = 1)),
      link = list(
        source = sources, target = targets,
        value =  as.numeric(confuse_matrix),
        color = colors
        )) %>% layout(title=session, font=list(family = "Arial",size = 20, color = 'black'))
}


confuse_bubblemat <- function(confuse_matrix_prop, label1, label2, session = "session")
{
prop <- as.numeric(confuse_matrix_prop)
data <- expand.grid(x = label1, y = label2) %>% bind_cols(prop = prop)
plot <- ggplot(data, aes(x = x, y = y, colour = prop, size = prop)) +
  geom_point()+
  scale_size_continuous(range = c(0, 10)) + 
  labs(x = "clusters", y = "inferred from") + theme_bw()

ggsave(paste0(session, ".svg"), plot = plot, device = svg, width = 5,height = 4)
}

## 返回最大的概率对应的index
func <- function(s, ident)
{
  if(max(s)>1.2/length(ident))
    return(ident[which(s == max(s))])
  else
    return("unassigned")
}

aero_colors_list <- as.character(lapply(colors_list, paste0, "A0")) #透明化颜色

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.